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ABSTRACT 


This  thesis  analyzes  the  dynamic  response  of  a  third  order  regulator  sys- 
tem. Particular  emphasis  is  placed  upon  the  loss  of  stability  of  the  nominal 
equilibrium  state.  The  system  utilized  in  this  research  models  the  funda- 
mental turning  dynamics  of  an  autonomous  vehicle.  We  make  extensive  use 
of  bifurcation  theory  methods  in  analyzing  the  dynamics  after  initial  loss 
of  stability.  The  effective  gain  of  the  system  is  used  as  the  main  bifurcation 
parameter,  since  this  is  directly  related  to  the  gain  margin  for  linear  systems. 
It  is  shown  that  the  nonlinear  characteristics  of  the  system  may  significantly 
affect  the  practical  significance  of  its  gain  margin,  as  a  measure  of  robustness 
to  parameter  variations,  unmodeled  dynamics,  and  external  disturbances. 
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I.  INTRODUCTION 


A.  PROBLEM  STATEMENT 

It  is  well  known  that  in  linear  dynamical  systems  one  of  the  most  popular 
ways  of  assessing  the  stability  properties  of  the  system  is  through  its  gain 
margin  (Friedland,  1986).  Roughly  speaking,  the  gain  margin  designates  the 
extent  at  which  the  effective  gain  of  the  system  can  be  increased  before  loss  of 
stability  occurs.  Therefore,  it  is  widely  used  in  linear  control  system  design 
in  order  to  quantify  a  measure  of  robustness  of  the  system  with  regards  to 
parameter  variations,  disturbances,  and  unmodelled  dynamics. 

In  this  work  we  examine  the  concept  of  gain  margin  in  the  light  of  nonlin- 
ear systems.  We  assume  that  the  baseline  linear  system  is  an  approximation 
to  a  nonlinear  system.  For  demonstration  purposes  we  employ  a  third  order, 
single  input  single  output  system,  with  cubic  nonlinearities.  This  system 
models  the  fundamental  turning  dynamics  of  a  marine  vehicle  (Oral,  1993). 
Primary  loss  of  stability  is  shown  to  occur  in  the  form  of  generic  bifurca- 
tions to  periodic  solutions  (Guckenheimer  k.  Holmes,  1983).  We  use  center 
manifold  reduction  techniques  and  integral  averaging  in  order  to  capture  the 
stability  properties  of  the  resulting  limit  cycles  (Chow  Sz  Mallet-Paret ,  1977). 
The  main  conclusion  of  this  work  is  that  the  linear  concept  of  a  gain  margin 
can  be  used  as  a  reliable  measure  of  stability  only  in  the  case  of  supercritical 
bifurcations  to  periodic  solutions.    In  the  case  of  subcritical  bifurcations,  a 


modification  is  necessary  which  is  based  on  the  nonlinear  characteristics  of 
the  system.  We  propose  the  use  of  a  parameter  which  governs  transitions 
from  supercritical  to  subcritical  bifurcations  as  a  nonlinear  gain  margin  of 
the  system.  Results  based  on  numerical  simulations  support  the  analytical 
predictions  of  this  work. 


B.  MATHEMATICAL  MODEL 

Consider  the  ideal  3rd-order  regulator  system  shown  in  Figure  1,  where 
y  is  the  actual  output,  j/ref  =  0  is  the  reference  input,  and, 

G(»)  =  "5 ^ .  (!) 

where  the  coefficients  a,-  correspond  to  a  stable  polynomial.   In  state  space 
form,  the  system  depicted  in  Figure  1  is  written  as, 

£,     =    x7  ,  (2) 

x7    =    x3  ,  (3) 

i3     =    -a7x3  -  aiir2  -  (a0  +  K)xx  ,  (4) 

where  the  state  vector  is, 

*i  =  y  i    «2  =  y  i    «3  =  y  •  (5) 

Physically,  we  can  think  of  this  system  as  a  representation  of  the  funda- 
mental turning  dynamics  of  a  marine  vehicle.  In  this  context,  Xy  represents 
the  vehicle's  lateral  deviation  from  the  commanded  straight  line  path,  x7  is 
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Figure  1:  Ideal  3rd-order  regulator  system 

the  orientation  angle,  and  x3  the  turning  rate.  To  account  for  the  geometric 
nonlinearities,  the  first  state  equation  (2)  is  written  as, 


xl  =x7  +  jx7  , 


(6) 


where  7  <  0  for  softening  spring  characteristics.  To  account  for  possible 
"over-steering"  or  "under-steering"  effects,  we  modify  the  output  of  the 
block  K  of  Figure  1,  to 


control    effort  =  —  Ky  —  K$yz 


(7) 


instead  of  —  Ky  of  the  linear  element.  Therefore,  the  nonlinear  system  under 

consideration  is, 

0  1  0         [  zi  1  -yx\ 

(8) 
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C.  OUTLINE  OF  ANALYSIS 

The  following  steps  are  performed  in  order  to  analyze  system  (8).  First, 
application  of  Routh's  criterion  yields  the  value  of  the  gain  margin,  or  the 
critical  value  of  K  for  stability.  Then  we  rewrite  the  system  of  equations  in  its 
normal  coordinate  form,  and  use  the  center  manifold  theorem  to  reduce  the 
system  into  a  two  dimensional  system.  We  apply  the  method  of  averaging 
to  the  reduced  system,  and  finally,  we  introduce  polar  coordinates  to  the 
averaged  system  in  order  to  reveal  the  existence  of  limit  cycles.  Development 
of  these  computations  is  the  subject  of  the  next  chapter. 


H.  ANALYSIS 

A.  LOSS  OF  STABILITY 

The  characteristic  equation  of  (8)  is, 

s3  +  a2s2  +a1s  +  (a0+K)  =0  .  (9) 

Application  of  Routh's  criterion  to  system  (9)  yields  the  critical  value  of  K, 
ifc,  for  stability  of  Xi  =  z7  =  x3  =  0, 

Kc  =  0=^2-00.  (10) 

If  K  <  Kc  the  system  is  stable,  whereas  for  K  >  Kc  it  becomes  unstable. 
The  value  of  KCi  given  in  equation  (10),  expressed  in  decibels  represents  the 
gain  margin  of  the  system. 

At  the  critical  point,  K  —  a\a2  —  a©,  the  characteristic  equation  (9)  is, 

sz  +  a2s7  +  Otis  +  aioza  =  0  , 

which  has  roots, 

*a  =  +\/5r« ,  (ii) 

Therefore,  we  can  see  that  the  above  loss  of  stability  is  characterized  by  the 
existence  of  a  pair  of  purely  imaginary  roots.  As  K  crosses  Kc,  one  pair  of 
complex  conjugate  roots  of  (9)  crosses  transversally  the  imaginary  axis. 
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A  situation  like  this  in  which  a  certain  parameter  is  varied  such  that 
the  real  part  of  one  pair  of  complex  conjugate  eigenvalues  of  the  linearized 
system  matrix  crosses  zero,  results  in  the  system  leaving  its  steady  state 
in  an  oscillatory  manner.  This  loss  of  stability  is  called  Hopf  bifurcation 
(Guckenheimer  &  Holmes,  1983)  and  generically  occurs  in  one  of  two  ways, 
supercritical  or  subcritical.  In  the  supercritical  case,  stable  limit  cycles  are 
generated  after  the  nominal  straight  line  motion  loses  its  stability.  The  am- 
plitudes of  these  limit  cycles  are  continuously  increasing  as  the  parameter 
distance  from  its  critical  value  is  increased.  For  small  values  of  this  critical- 
ity  distance  the  resulting  limit  cycle  is  of  small  amplitude  and  differs  little 
from  the  initial  nominal  state.  In  the  subcritical  case,  however,  stable  limit 
cycles  are  generated  before  the  nominal  state  loses  its  stability.  Therefore, 
depending  on  the  initial  conditions  it  is  possible  to  diverge  away  from  the 
nominal  straight  line  path  and  converge  towards  a  limit  cycle  even  before  the 
nominal  motion  loses  its  stability.  This  means  that  in  the  subcritical  Hopf 
bifurcation  case  the  domain  of  attraction  of  the  nominal  state  is  decreasing 
and  in  fact  it  shrinks  to  zero  as  the  critical  point  is  approached.  Random 
external  disturbances  of  sufficient  magnitude  can  throw  the  vehicle  off  to  an 
oscillatory  steady  state  even  though  the  nominal  state  may  still  remain  sta- 
ble. After  the  nominal  state  becomes  unstable,  a  discontinuous  increase  in 
the  magnitude  of  motions  is  observed  as  there  exist  no  simple  stable  nearby 
attractors  for  the  vehicle  trajectory  to  converge  to.  Distinction  between  these 
two  qualitatively  different  types  of  bifurcation  is,  therefore,  essential  in  the 


design  of  the  autopilot.  The  computational  procedure  requires  examination 
of  the  higher  order  terms  in  the  equations  of  motion  and  it  is  the  subject  of 
the  next  section. 


B.  COORDINATE  TRANSFORMATIONS 

System  (8)  is  written  in  the  form, 

x  =  Ax  +  g(x)  , 


(12) 


where  A  is  the  linearized  system  matrix  and  g(x)  contains  the  cubic  terms. 
At  the  bifurcation  point  K  =  Kc,  matrix  A  has  eigenvalues  given  by  (11), 
and  eigenvectors, 


0  l/v^T     l 

1  0         -a2 

0     -y/5[      a\ 


(13) 


By  taking  T  to  be  the  matrix  of  critical  eigenvectors  of  A,  the  transformation, 


x  =  Tz  ,     z  =  T_1x  , 

transforms  system  (12)  into  its  normal  coordinate  form, 
z  =  T"1ATz+T-1g(Tz), 


where, 


and, 


T1 


<*!  +a| 


T_1AT  = 


«1  «2         <*1  +  <*2  "2 


■y/*l 
1 


0       -yjal      0 
oTT        0  0 

0  0         -a2 


(14) 


(15) 


(16) 


(17) 


For"  values  of  K  close  to  the  bifurcation  point  Ke,  the  matrix  T  *AT 
becomes, 


T-iAT 


a'e  -{y/al  +  uJe)  0 

^/oT  +  Je  a'e  0 

0  0  -a2  +  p'e 


(19) 


where  e  denotes  the  criticality  difference, 

e  =  K-Kc, 

and  a'  is  the  derivative  of  the  real  part  of  the  critical  eigenvalues  with  respect 
to  e,  u;'  is  the  derivative  of  the  imaginary  part  of  the  critical  eigenvalues  with 
respect  to  e,  and  j/  is  the  derivative  of  the  third  eigenvalue  53  with  respect 
to  e.  These  are  computed  using  a  perturbation  series  approach  as  follows. 
The  characteristic  equation  (9)  can  be  written  in  the  form, 

a3  +  a2*2+a1«  +  (a1a2+e)=0,  (20) 

where  we  have  used  (10)  and  (19).  The  roots  of  (20)  are  expressed  as, 

3X    =    a'e-(^  +  u'e)i1  (21) 

32     =    a'e  +  {^  +  u>'e)i  ,  (22) 

s3    =    -a2  +  p'e  •  (23) 

K  we  substitute  (21)  through  (23)  in  (20)  and  neglect  terms  of  order  e2,  e3, 
we  get, 

(24) 

(25) 

(26) 


2(ai+a!)' 

2v^T(a1+a22)  ' 
1 


P     = 


8 


C.  REDUCTION  OF  ORDER 

The  physical  variables  x,  are  related  to  the  normal  coordinates  z±  through 
(14),  and  using  (13)  we  get, 

1 

x2     =    zx  —  a2Z3  ,  (27) 

a*     =     -y/o^zi  +  «|^3  • 

It  can  be  seen  from  (26)  that  j/  <  0,  and  therefore  the  eigenvalue  j3  is 
locally  (for  e  small)  negative,  as  shown  in  equation  (23).  In  fact,  a  quick 
root  locus  plot  of  (9)  will  show  that  33  is  negative  for  all  values  of  K  > 
0.  Therefore,  the  flow  of  (8)  in  the  direction  of  z$  converges  to  zero.  All 
interesting  bifurcation  phenomena  are  locally  restricted  on  a  two  dimensional 
manifold  that  decsribes  the  time  evolution  of  the  critical  coordinates  Z\ ,  z2 ; 
this  is  the  center  manifold  of  (8).  According  to  the  center  manifold  theorem 
(Guckenheimer  &  Holmes,  1983),  the  stable  coordinate  z$  can  be  expressed 
as  a  function  of  the  critical  coordinates  z1 ,  z7 ,  and  this  relationship  is  at  least 
of  quadratic  order.  In  fact  due  to  the  symmetry  in  our  problem,  the  above 
relationship  is  of  third  order.  Therefore,  23  does  not  affect  the  nonlinear 
terms  in  (15)  and  we  can  write  (27)  in  the  form, 

1 

x7    =    zx  ,  (28) 

x3     =     -\/ai>2  • 


Using  (28)  we  can  write, 

g(Tz) 
and, 


iA 


3/2 


T-^Tz) 


«i  +ai 


I  -K34/a\ 


a1a2yzf  -  K3z%a2/ai 
al7zf  -  K*z\l<xT 


(29) 


Using  equations  (18)  and  (29)  we  substitute  in  (15)  and  write  the  normal 
equations  in  Zi,  z2  as, 


ii     =    a'ez!  -  {y/al  +  dz)z2  +  F1(zliz2)  , 
z2     =    (x/aT  +  u>'e:)zi  +  o/ez2  +  ^2(21,22)  > 


where, 


r/  \  aia2         3 


"2 


(a1+oi)a?/2 


K$z^  , 


K*z« 


(30) 
(31) 


(32) 
(33) 


a2  +0$  '   *       (ax  +a%)a1^3  2  * 
Equations  (30)  and  (31)  describe  the  suspended  center  manifold  flow  of  (15) 
(Guckenheimer  &  Holmes,  1983). 


D.  AVERAGING 

If  we  introduce  polar  coordinates  in  the  form, 

z1  =  R  cos  9  ,      Z2  =  R  sin  0  , 


(34) 
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we  can  write  equations  (30),  (31)  as, 

R    =    a' eR  +  FX(R  cos  B,R  sin  B)  cos  B 

+F2  (R  cos  0,  R  sin  B)  sin  9  ,  (35) 

R9    =    (y/a~^  +  u' e) R +  F2(R  cos  B,R  an  B)  cos  B 

-Fx  (R  cos  9,  R  sin  0)  sin  9  .  (36) 

Equation  (35)  is  written  in  the  form, 

R  =  a'eR  +  V{B)&  ,  (37) 

where  the  function  V{9)  is  27r-periodic  in  the  angular  coordinate  9, 

V{9)    =       aiQ!227cos4fl "2      3/7K3cosBsiii3B 


V^l 34,..    ,     1 

ax(ax  +a\) 


7<ns39nn9+—T—r-j^K3&nA9.      (38) 


If  equation  (37)  is  averaged  over  one  cycle  in  0,  we  get  an  equation  with 
constant  coefficients, 

R  =  a'eR  +  QR3  ,  (39) 

where 

Q"iC^')dt-  (40) 

Substituting  (38)  into  (40)  and  evaluating  the  integral  yields, 


3    ala7~r  +  K3 
8'  aii^+al) 


q=^?z':z-  (4i) 


Similar  averaging  can  be  performed  for  equation  (36)  which  has  the  form, 
B  =  V^T  +  u/e  +  T{B)R2  ,  (42) 
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where    - 

T(6)     =      y^a\1cos/ie+——^ rrtfacosflsin3* 

w  a1+a22I  ^(ai+ai) 

^V/  cos3  6  sin  0+     ...    a2 K3anU,     (43) 

and  we  have  assumed  R  ^  0.  The  averaged  form  of  (42)  is, 

6  =  yfc  +  u/e  +  GR7  ,  (44) 


where 


and  using  (43), 


"■sf'w*'  <«> 


g  =  ^-  X'     .-■  («) 


3a2      aja27  —  jFT3 
8    '«?/'(«»,+ a?) 

The  system  of  equations  (39)  and  (44)  exhibits  similar  stability  properties 

to  the  original  system  (37)  and  (42)  (Chow  &  Mallet-Paret,  1977),  and  is 
studied  in  the  following  section. 


E.  LIMIT  CYCLE  ANALYSIS 

Equation  (39)  has  two  steady  state  solutions,  one  at  R  =  0  which  corre- 
sponds to  the  trivial  equilibrium  solution  at  zero,  and  one  at 

*S  -  ~*  •  (47) 

This  equilibrium  solution  corresponds  to  a  periodic  solution  or  limit  cycle 
in  the  cartesian  coordinates  z1?  z-i  from  (64).    Since  a'  as  seen  by  (24)  is 
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always,  positive,  existence  of  these  periodic  solutions  depends  on  the  value  of 
Q.  Specifically, 

•  if  Q  <  0,  periodic  solutions  exist  for  e  >  0  or  K  >  Kc,  and 

•  if  Q  >  0,  periodic  solutions  exist  for  e  <  0  or  K  <  Kc. 

The  Floquet  exponent  of  (39)  in  the  vicinity  of  (47)  is 

P  =  -2a'e  ,  (48) 

and  we  can  see  that 

•  if  periodic  solutions  exist  for  K  >  Kc  they  are  stable,  and 

•  if  periodic  solutions  exist  for  K  <  Kc  they  are  unstable. 

We  refer  to  the  first  case  as  the  supercritical  Poincare-Andronov-Hopf  (PAH) 
bifurcation  and  to  the  second  case  as  the  subcritical  PAH  bifurcation  (Guck- 
enheimer  &  Holmes,  1983). 


The  period  of  these  limit  cycles  is  computed  by  substituting  (47)  in  (44), 

T=  ^w£+gflr^(1-^rH+0(e)-    (49) 

The  amplitude  of  the  limit  cycles  is  computed  from  (47)  and  (28),  and  in 
terms  of  the  physical  variables  Xi,  x^,  jt3,  we  get, 


_       \/3    laia2  -  op  ^K 
Xl     ~        3  V     «la27  +  K3     ' 


(50) 


ttiCojOj  -  oo  -  K) 


3  \        a{<x2i  +  Kz 
13 


We  can  see  that  in  the  supercritical  case,  on  loss  of  stability  of  equilibrium 
the  steady  state  becomes  a  periodic  oscillatory  state,  the  amplitude  of  the 
oscillation  being  proportional  to  the  square  root  of  the  criticality  e,  the  differ- 
ence of  the  gain  K  from  its  critical  value  Kc  at  which  stability  of  equilibrium 
is  lost.  This  form  of  loss  of  stability  is  called  "soft"  loss  of  stability  since 
the  oscillating  state  for  small  e  differs  little  from  the  equilibrium  state.  In 
the  subcritical  case,  before  the  steady  state  loses  stability  the  domain  of 
attraction  becomes  very  small  as  is  bounded  by  the  amplitudes  of  the  un- 
stable limit  cycles,  and  a  random  disturbance  can  throw  the  system  off  its 
equilibrium  state  even  before  its  domain  of  attraction  has  completely  disap- 
peared. This  form  of  loss  of  stability  is  called  "hard" .  Here  the  system  leaves 
its  steady  state  with  a  jump  to  a  different  state  of  motion  which  can  be  a 
stable  oscillation  with  a  locally  discontinuous  increase  in  the  amplitude,  a 
more  complicated  bounded  motion,  or  even  an  unbounded  motion  depending 
on  other  higher-than-third  order  terms  that  have  not  been  incorporated  in 
system  (8). 
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m.  RESULTS  AND  DISCUSSION 


A.  RESULTS 

For  demonstration  purposes  we  assume  that  in  the  case  K  =  0  the  system 
has  been  designed  in  accordance  to  the  ITAE  optimal  criterion  for  third  order 
systems;  i.e., 

a2     =     1.75a;n 

«!    =    2.15u£,  (53) 

«o    =    «J , 

where  the  natural  frequency  u;n  is  in  general  selected  according  to  the  desired 
bandwidth  of  the  system.  The  ITAE  criterion  is  a  standard  performance 
index  and  it  minimizes  the  integral  time  absolute  error  /0°°  t\e(t)\  dt,  of  the 
step  response  e(t)  of  a  system  (Dorf,  1992).  This  criterion  is  satisfied  for 
a  third  order  system  when  the  coefficients  of  its  characteristic  equation  are 
selected  as  in  equations  (53).  The  critical  value  of  K  for  stability  is  computed 

by, 

Kc  =  axa7  -  oo  =  2.7625u;J  .  (54) 

A  graphical  representation  of  equation  (54)  is  shown  in  Figure  2  for  a  nominal 
range  of  u)n  between  1  and  2  rad/sec. 

The  cubic  coefficient  Q  that  dictates  the  nature  (supercritical  or  subcrit- 
ical)  and  stability  of  the  resulting  limit  cycles  as  K  exceeds  Kc  is  given  by 
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Figure  2:  Critical  value  Kc  versus  un 
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Figure  3:  Cubic  coefficient  Q  versus  cvn  for  A'3  =  -5  and  different  values  of  7 
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(41),  and  using  (53)  we  get, 

Q  =  0.2707u>n7  +  0.0112-f-  .  (55) 

A  plot  of  equation  (55)  for  K3  =  —5  and  four  different  values  of  7  is  shown 
in  Figure  3. 

Based  on  the  results  presented  in  Figure  3,  the  following  observations  can 
be  made:  (1)  For  fixed  Ks  and  a  given  value  of  u>n,  supercritical  bifurca- 
tions are  ensured  for  a  value  of  7  less  than  a  critical  threshold,  computed  by 
equating  Q  =  0.  (2)  In  the  case  of  supercritical  bifurcations,  we  expect  to  see 
an  oscillatory  response  approaching  zero  when  K  <  Kc  and  converging  to  a 
periodic  solution  for  K  >  Kc.  The  above  response  should  be  independent  of 
the  initial  conditions,  at  least  locally.  (3)  In  the  case  of  subcritical  bifurca- 
tions, we  expect  to  see  an  oscillatory  response  for  K  <  Kc  which  may  or  may 
not  converge  to  zero.  This  depends  on  the  initial  conditions.  If  the  response 
diverges  from  zero,  the  final  attractor  can  be  another  oscillation  or  a  more 
complicated  motion.  The  same  is  true  for  K  >  Ke,  here  the  response  should 
diverge  regardless  of  the  initial  conditions.  These  conclusions  are  confirmed 
in  the  next  section  using  direct  numerical  integrations  of  equations  (8). 


B.  NUMERICAL  SIMULATIONS 

An  example  of  supercritical  behavior  is  shown  in  Figures  4  and  5.  In  both 
of  these  figures  we  present  results  based  on  direct  numerical  integrations  of 


17 


equations  (8)  for  the  following  conditions, 

wn  =  2  ,     K  =  1.05^  ,     K3  =  -5  ,     7  =  -0.1  .  (56) 

For  these  conditions,  which  are  selected  to  demonstrate  the  supercritical  case, 
we  can  see  from  Figure  3  that  Q  <  0  and,  therefore,  we  have  a  supercritical 
Hopf  bifurcation.  Figure  4  shows  clearly  the  development  of  a  stable  periodic 
solution  with  amplitude  x7  =  0.68  and  period  T  =  2.14.  The  theoretical 
value  for  the  limit  cycle  amplitude  is  computed  from  (56)  and  (51),  as  x2  — 
0.6405  which  is  in  excellent  agreement  with  the  numerical  value.  Likewise, 
the  period  of  oscillation  is  found  from  (56)  and  (49),  as  T  =  2.1351  which  is 
also  close  to  the  actual  period.  Figure  5  shows  the  convergence  to  the  limit 
cycle  in  the  (3:25*2)  phase  subspace  using  two  sets  of  initial  conditions.  One 
set  (xx,x2,X3)  =  (0,0.5,0)  is  located  inside  the  limit  cycle  and  the  other  set 
(xi,x2,x3)  =  (0,1,0)  is  located  outside  the  limit  cycle.  It  can  be  seen  that 
both  trajectories  converge  to  the,  numerically  computed,  periodic  solution. 

An  example  of  subcritical  behavior  for  K  >  K c  is  shown  in  Figures  6  and 
7.  In  both  of  these  figures  the  conditions  were  the  same  as  (56)  with  the 
exception  of  7  which  was  7  =  0.2  for  Figure  6  and  7  =  0.1  for  Figure  7. 
These  correspond  to  Q  >  0  as  can  be  seen  from  Figure  3,  which  confirms  the 
subcritical  behavior.  The  motion  which  corresponds  to  7  =  0.2  is  strongly 
subcritical  and  it  becomes  unbounded  shortly  after  20  seconds.  The  motion 
which  corresponds  to  7  =  0.1  is  less  subcritical  since  it  yields  a  smaller,  but 
still  positive,  value  for  Q.   As  seen  in  Figure  7,  it  appears  that  trajectories 
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Figure  4:  Time  history  {t,x2)  for  ujn  =  2,  K  =  1.05J3fe,  #3  =  -5,  and  7  =  -0.1 


Figure  5:  Phase  subspace  plot  (x2,X2)  for  wn  =  2,  K  =  1.05KC,  K3  t=  -5,  7  =  -0.1, 
and  two  initial  conditions 
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Figure  6:  Time  history  {t,z2)  for  wn  =  2,  K  =  1.05ATC,  Kz  =  -5,  and  j  =  0.2 
converge  to  a  relaxation  oscillation  of  high  amplitude. 

An  example  of  subcritical  behavior  for  K  <  Kc  is  shown  in  Figures  8  and 
9,  for  the  following  conditions, 


u>n  =  2  ,     K  =  0.95KC  ,     K3  =  -5  ,     7  =  0.2  , 


(57) 


and  different  initial  conditions  in  x2.  For  small  initial  conditions,  (x1,x2,x3)  = 
(0,0.1,0),  the  system  is  located  inside  the  unstable  limit  cycle  and  it  con- 
verges to  the  stable  equilibrium,  as  shown  in  Figure  8.  However,  for  large 
initial  conditions,  (xx,x7,X3)  =  (0,0.5,0),  the  system  trajectory  is  located 
outside  the  unstable  limit  cycle  and  it  diverges  away  from  the  equilibrium, 
even  though  this  equilibrium  is  still  stable.  It  appears  from  Figure  9  that,  for 
small  7,  the  system  converges  to  a  relaxation  oscillation,  while  higher  value 
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Figure  7:  Phase  subspace  plot  (x2,*2)  for  u)n  —  %  K  =  1.05iifc,  #3  =  -5,  and 
7  =  0.1 

of  7  generate  more  severe  subcritical  behavior  and  the  motion  becomes  un- 
bounded, as  before.  In  the  case  of  supercritical  behavior  and  for  K  <  Kc, 
we  observed  that  numerical  integrations  converged  to  the  stable  equilibrium 
regardless  of  the  initial  condition  in  x2,  as  they  should. 


C.  MULTIPLE  EQUILIBRIUM  STATES 

So  far,  our  analysis  has  been  on  stability  properties  of  the  trivial  equilib- 
rium of  (8),  and  its  bifurcations  to  periodic  solutions.  It  is  possible,  however, 
that  additional  equilibrium  points  may  exist.   To  explore  this,  we  write  sys- 
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Figure  8:    Time  history  {t,x2)  for  wn  =  2,  K  =  0.95A"C,  K3  = 
x2(t  =  0)  =  0.1 


-5,  7  =  0.2,  and 


Figure  9:    Time  history  (t,x7)  for  con  =  2,  K  -  0.95 Kc,  K3  =  -5,  7  =  0.2,  and 
x2(t  =  0)  =  0.5 
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Figure  10:  Steady  state  solutions  X\  versus  K3  for  K  —  0.95iTc,  7  =  0,  and  three 
values  of  u>„ 


tern  (8)  in  the  form, 


xi     =    x7  +  7^2  « 

x2     =    x3  , 

x3     =     -(a0  +  K)xx  -  axx2  -  oc7x3  -  Kzx\ 


(58) 
(59) 
(60) 


In  order  to  compute  its  equilibrium  points  we  must  set  the  time  derivatives 
ii,  ^2,  &3  zero,  and  solve  for  the  equilibrium  point  Xi,  X7,  X3.  We  examine 
the  case  7  >  0  first. 

Equation  (58)  yields, 

X2+7X23=0,  (61) 

which,  since  7  >  0,  has  only  one  solution,  namely  X2  =  0.    Equation  (59) 
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Figure  11:  Time  history  (t,xi)  for  u>n  =  2,  K  =  0.95Jifc,  iJT3  =  -5,  7  =  0,  and  two 
initial  conditions  in  x\ 

gives  X3  =  0,  and  then  (60)  can  be  solved  for  the  remaining  equilibrium 
solution  Xi, 

Xtfao  +  K  +  KZX\)  =  0  .  (62) 

Elquation  (62)  admits  the  trivial  solution,  X\   =  0,  always,  and  two  more 
possible  solutions  provided  K3  <  0,  given  by, 

(63) 


K3 

Equation  (63)  yields  two  additional  symmetrically  located  equilibrium  points 
which  are  generated  as  K3  becomes  negative.  A  typical  plot  for  K  =  0.95KC 
and  three  values  of  a>n  =  1,  1.5,  2  is  presented  in  Figure  10.  Therefore,  we  see 
that  the  same  coefficient,  if3,  which  governs  the  transition  from  subcritical 
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to  supercritical  bifurcations  to  periodic  solutions  of  the  trivial  equilibrium 
state,  is  also  associated  with  the  existence  of  additional  equilibrium  states. 

In  order  to  analyze  the  stability  properties  of  these  solutions  we  linearize 
equations  (58)  through  (60)  in  the  neighborhood  of  X;.  If  we  denote  by  &  the 
deviation  of  Xi  from  equilibrium;  i.e,  £  =  sc,-  —  X;,  we  can  write  the  linearized 
system  as, 

6=6, 

6    =    6,  (64) 

6    =    -(ao  +  i06-«i6-«26-  3*3^6  • 

If  we  substitute  (63)  in  (64)  we  get, 

6=6, 

6=6,  (65) 

6     =    2(a0  +  ii:)6- ^6-026. 

The  characteristic  equation  of  (65)  is, 

s3 +a2s2 -\-a1s-2(a0+K)  =  0,  (66) 

which  means  that  the  additional  non-trivial  equilibria  are  clearly  unstable 
with  divergent  dynamics.  Therefore,  the  case  K3  <  0  which  was  shown  to 
be  beneficial  from  the  point  of  view  of  Hopf  bifurcations,  is  undesirable  from 
the  point  of  view  of  static  bifurcations.  The  final  compromise  depends  of 
course  on  the  particular  demands  and  specifications  of  the  design.  Figure  11 
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Figure  12:  Steady  state  solutions  X\  versus  K3  for  K  —  d.9bKe,  7  =  -0.1,  u>n  =  2, 
and  X7  ^  0 

demonstrates  the  divergent  properties  of  the  additional  equilibrium  states 
for  K3  <  0.  Two  numerical  simulations  (t,Xi)  are  shown  for  Kz  =  —5, 
ion  =  2,  and  if  =  0.95iiTc>  and  for  two  sets  of  initial  conditions,  (ei,^?^)  = 
(2.40,0,0)  and  (xi,X2,x3)  =  (2.41,0,0).  The  unstable  equilibrium  Xi  is 
located  at  2.408  as  predicted  by  (63).  We  can  see  that,  as  expected,  numerical 
simulations  inside  the  stable  potential  well,  xx  <  X\ ,  converge  to  zero  while 
those  outside,  Xi  >  X\t  diverge  and  quickly  become  unbounded. 

The  case  7  <  0  can  be  analyzed  similarly.   Equation  (61)  admits  in  this 
case  the  trivial  solution  X7  =  0  which  was  analyzed  previously,  and  two  more 
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Figure  13:   Degree  of  stability  of  the  lower  outer  solutions  which  appear  in  Figure 
12 


Figure  14:  Trajectories  for  K3  =  -4.6  and  two  different  initial  conditions;  variables 
correspond  to  Figure  12 
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Figure  15:    Time  history  (t,r,i)  for  A3  =  -3  and  two  initial  conditions;  variables 
correspond  to  Figure  12 
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symmetrically  located  solutions, 

XI  =  -i  .  (67) 

Then  X3  =  0  while  Xi  is  obtained  from  (60)  as  the  solution  to, 

(a0  +  K)Xt  +  axX2  +  KZX\  =  0  .  (68) 

A  typical  solution  set  of  (68)  is  presented  in  Figure  12  for  u;n  =  2,  K  = 
0.95Ke,  and  7  =  —0.1  with  X7  given  from  (67).  If  we  compare  this  to  the 
corresponding  solutions  for  7  =  0  shown  in  Figure  10  we  can  observe  that 
the  nontrivial  solutions  for  7  =  0  retain  their  shape,  these  are  denoted  by 
outer  in  Figure  12.  The  main  difference  here  is  the  existence  of  two  outer 
solutions  for  most  of  the  range  of  K3.  The  trivial  solution  Xj  =  0  of  Figure 
10  perturbs  in  this  case  into  the  solution  labeled  as  inner  in  Figure  12.  It 
can  also  be  observed  that  there  is  a  point,  in  this  case  at  about  K3  =  —4.8 
where  the  lower  outer  solution  coalesces  with  the  inner  solution.  For  values 
of  K3  less  than  this  critical  point  only  the  upper  outer  solution  remains. 

Stability  properties  of  these  solutions  can  be  established  by  linearization. 
The  linearized  system  of  (58)  through  (60)  for  X7  /  0  takes  the  form, 

6    =   -26, 

6=6,  (69) 

6     =    -(a0  +  K  +  ZK3X*)t1-a1Z2-a2k. 

The  characteristic  equation  of  (69)  is, 

s*  +  a2s7  +  alS  -  2(a0  +  K  +  3K3X7)  .  (70) 
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Numerical  computation  of  the  roots  of  (70)  revealed  that  both  the  inner  solu- 
tion and  the  upper  outer  solutions  were  unstable.  The  lower  inner  solution  is 
initially  stable,  in  other  words  the  point  where  the  lower  outer  and  the  inner 
solutions  of  Figure  12  meet  is  a  saddle-node  point.  The  degree  of  stability 
of  the  lower  outer  solution,  defined  as  the  largest  real  part  of  the  three  roots 
of  (70),  is  plotted  in  Figure  13  versus  K3.  This  corresponds  to  a  complex 
conjugate  root.  We  can  see  that  it  is  initially  stable  and  that  it  undergoes  a 
Hopf  bifurcation  at  a  value  of  K3  approximately  —4.15. 

These  results  are  confirmed  by  the  numerical  integrations  presented  in 
Figures  14  through  16.  For  K3  =  —4.6,  the  lower  outer  equilibrium  solution 
is  stable.  Therefore,  the  numerical  simulations  will  converge  to  either  this 
solution  or  the  trivial  equilibrium  depending  on  the  initial  conditions,  see 
Figure  14.  The  trivial  equilibrium  is  a  rather  weak  attractor  in  this  case 
since  K  is  very  near  its  critical  value  Kc.  For  Kz  =  —3,  the  lower  outer 
equilibrium  point  has  become  unstable  and  is  surrounded  by  a  stable  Hmit 
cycle.  Depending  on  the  initial  conditions,  trajectories  will  converge  to  either 
the  stable  trivial  equilibrium  or  one  of  the  two  stable  limit  cycles,  see  Figure 
15.  Apparently  there  exists  a  region  where  the  final  outcome  is  sensitive 
to  the  choice  of  initial  conditions,  and  the  corresponding  transient  response 
resembles  a  random  pre— chaotic  motion  before  the  trajectories  converge  to 
the  corresponding  attractor.  This  limit  cycle  persists  as  K  exceeds  Kc,  see 
Figure  16.  In  this  case  the  trajectories  will  converge  to  either  the  trivial  or 
the  non-trivial  limit  cycles,  depending  again  on  the  initial  conditions. 
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IV.  CONCLUSIONS  AND  RECOMMENDATIONS 


This  work  presented  a  methodology  for  assessing  the  dynamic  response  of 
a  third  order  system  with  respect  to  changes  in  its  gain.  Choice  for  the 
baseline  system  was  motivated  by  modeling  the  fundamental  dynamics  of  an 
autonomous  vehicle.  The  methods,  however,  are  of  general  nature  and  can 
be  applied  to  any  given  system.  An  extensive  study  of  the  dynamic  loss  of 
stability  was  performed  based  on  Hopf  bifurcation  theory  techniques.  The 
existence  of  both  subcritical  and  supercritical  bifurcations  to  periodic  solu- 
tions was  established  depending  on  the  system  parameters.  It  was  shown  that 
the  critical  system  gain  for  stability  is  useful  for  design  purposes  only  to  the 
extent  that  it  is  accompanied  by  supercritical  bifurcations.  The  latter  can  be 
studied  using  a  comprehensive  nonlinear  study  like  the  one  presented  in  this 
work.  Future  work  should  concentrate  on  classifying  the  nonlinear  dynamics 
of  various  dynamical  systems  in  terms  of  their  order  and  linear/nonlinear 
properties.  Such  a  generalization  could  allow  the  establishment  of  a  more 
reliable  set  of  measures  than  the  linear  gain  margin,  which  would  be  used  in 
control  system  design. 
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